# Calculate the typical values for all variables in Model 3, except for quota,
# which we'll force to be FALSE and TRUE
data_to_predict <- typical(model = model_frac_logit3,
quota = c(FALSE, TRUE))
ggplot(frac_logit_pred, aes(x = quota, y = predicted)) +
geom_pointrange(aes(ymin = conf.low, ymax = conf.high)) +
labs(x = "Quota", y = "Predicted proportion of women in parliament") +
theme_clean()
# Calculate the typical values for all variables in Model 3, except for quota,
# which we'll force to be FALSE and TRUE
data_to_predict <- data_grid(model = model_frac_logit3,
quota = c(FALSE, TRUE))
data_to_predict
# Calculate the typical values for all variables in Model 3, except for quota,
# which we'll force to be FALSE and TRUE
data_to_predict <- datagrid(model = model_frac_logit3,
quota = c(FALSE, TRUE))
data_to_predict
# Plug this new dataset into the model
frac_logit_pred <- predictions(model_frac_logit3, newdata = data_to_predict)
frac_logit_pred
ggplot(frac_logit_pred, aes(x = quota, y = predicted)) +
geom_pointrange(aes(ymin = conf.low, ymax = conf.high)) +
labs(x = "Quota", y = "Predicted proportion of women in parliament") +
theme_clean()
data_to_predict <- typical(model = model_frac_logit3,
polyarchy = seq(0, 100, by = 1))
data_to_predict <- datagrid(model = model_frac_logit3,
polyarchy = seq(0, 100, by = 1))
frac_logit_pred <- predictions(model_frac_logit3,
newdata = data_to_predict)
ggplot(frac_logit_pred, aes(x = polyarchy, y = predicted)) +
geom_ribbon(aes(ymin = conf.low, ymax = conf.high),
alpha = 0.4, fill = "#480B6A") +
geom_line(size = 1, color = "#480B6A") +
labs(x = "Polyarchy (democracy)",
y = "Predicted proportion of women in parliament") +
theme_clean()
summary(frac_logit_pred)
View(frac_logit_pred)
ggplot(frac_logit_pred, aes(x = polyarchy, y = estimate)) +
geom_ribbon(aes(ymin = conf.low, ymax = conf.high),
alpha = 0.4, fill = "#480B6A") +
geom_line(size = 1, color = "#480B6A") +
labs(x = "Polyarchy (democracy)",
y = "Predicted proportion of women in parliament") +
theme_clean()
ggplot(frac_logit_pred, aes(x = quota, y = estimate)) +
geom_pointrange(aes(ymin = conf.low, ymax = conf.high)) +
labs(x = "Quota", y = "Predicted proportion of women in parliament") +
theme_clean()
# Calculate the typical values for all variables in Model 3, except for quota,
# which we'll force to be FALSE and TRUE
data_to_predict <- datagrid(model = model_frac_logit3,
quota = c(FALSE, TRUE))
data_to_predict
# Plug this new dataset into the model
frac_logit_pred <- predictions(model_frac_logit3, newdata = data_to_predict)
frac_logit_pred
ggplot(frac_logit_pred, aes(x = quota, y = estimate)) +
geom_pointrange(aes(ymin = conf.low, ymax = conf.high)) +
labs(x = "Quota", y = "Predicted proportion of women in parliament") +
theme_clean()
marginaleffects(model_frac_logit3, variables = "polyarchy") %>%
tidy()
marginaleffects(model_frac_logit3,
newdata = typical(polyarchy = c(56, 66)))
modelsummary(list("OLS" = model_ols3,
"Fractional logit" = model_frac_logit3,
"Fractional logit<br>(marginal effects)" = marginaleffects(model_frac_logit3)),
gof_map = gof_stuff,
escape = FALSE)
model_beta <- betareg(prop_fem ~ quota | quota,
data = vdem_2015,
link = "logit")
vdem_2015_fake0 <- vdem_2015 %>%
mutate(prop_fem = ifelse(prop_fem == 0, 0.001, prop_fem))
model_beta <- betareg(prop_fem ~ quota | quota,
data = vdem_2015_fake0,
link = "logit")
tidy(model_beta)
beta_mu_intercept <- model_beta %>%
tidy() %>%
filter(component == "mean", term == "(Intercept)") %>%
pull(estimate)
beta_mu_quota <- model_beta %>%
tidy() %>%
filter(component == "mean", term == "quotaTRUE") %>%
pull(estimate)
plogis(beta_mu_intercept + beta_mu_quota) - plogis(beta_mu_intercept)
beta_phi_intercept <- model_beta %>%
tidy() %>%
filter(component == "precision", term == "(Intercept)") %>%
pull(estimate)
beta_phi_quota <- model_beta %>%
tidy() %>%
filter(component == "precision", term == "quotaTRUE") %>%
pull(estimate)
no_quota_title <- paste0("dprop(mean = plogis(", round(beta_mu_intercept, 2),
"), size = exp(", round(beta_phi_intercept, 2), "))")
quota_title <- paste0("dprop(mean = plogis(", round(beta_mu_intercept, 2),
" + ", round(beta_mu_quota, 2),
"), size = exp(", round(beta_phi_intercept, 2),
" + ", round(beta_phi_quota, 2), "))")
ggplot(data = tibble(x = 0:1), aes(x = x)) +
geom_density(data = vdem_2015_fake0,
aes(x = prop_fem, fill = quota), alpha = 0.5, color = NA) +
stat_function(fun = dprop, size = 1,
args = list(size = exp(beta_phi_intercept),
mean = plogis(beta_mu_intercept)),
aes(color = no_quota_title)) +
stat_function(fun = dprop, size = 1,
args = list(size = exp(beta_phi_intercept + beta_phi_quota),
mean = plogis(beta_mu_intercept + beta_mu_quota)),
aes(color = quota_title)) +
scale_fill_viridis_d(option = "plasma", end = 0.8, name = "Quota",
guide = guide_legend(ncol = 1, order = 1)) +
scale_color_viridis_d(option = "plasma", end = 0.8, direction = -1, name = "",
guide = guide_legend(reverse = TRUE, ncol = 1, order = 2)) +
labs(x = "Proportion of women in parliament") +
theme_clean() +
theme(legend.position = "bottom")
model_beta_no_phi <- betareg(prop_fem ~ quota,
data = vdem_2015_fake0,
link = "logit")
tidy(model_beta_no_phi)
ggplot(data = tibble(x = 0:1), aes(x = x)) +
geom_density(data = vdem_2015_fake0,
aes(x = prop_fem, fill = quota), alpha = 0.5, color = NA) +
stat_function(fun = dprop, size = 1,
args = list(size = beta_phi,
mean = plogis(beta_mu_intercept)),
aes(color = no_quota_title)) +
stat_function(fun = dprop, size = 1,
args = list(size = beta_phi,
mean = plogis(beta_mu_intercept + beta_mu_quota)),
aes(color = quota_title)) +
scale_fill_viridis_d(option = "plasma", end = 0.8, name = "Quota",
guide = guide_legend(ncol = 1, order = 1)) +
scale_color_viridis_d(option = "plasma", end = 0.8, direction = -1, name = "",
guide = guide_legend(reverse = TRUE, ncol = 1, order = 2)) +
labs(x = "Proportion of women in parliament") +
theme_clean() +
theme(legend.position = "bottom")
beta_phi <- model_beta_no_phi %>%
tidy() %>%
filter(component == "precision") %>%
pull(estimate)
no_quota_title <- paste0("dprop(mean = plogis(", round(beta_mu_intercept, 2),
"), size = ", round(beta_phi, 2), ")")
quota_title <- paste0("dprop(mean = plogis(", round(beta_mu_intercept, 2),
" + ", round(beta_mu_quota, 2),
"), size = ", round(beta_phi, 2), ")")
ggplot(data = tibble(x = 0:1), aes(x = x)) +
geom_density(data = vdem_2015_fake0,
aes(x = prop_fem, fill = quota), alpha = 0.5, color = NA) +
stat_function(fun = dprop, size = 1,
args = list(size = beta_phi,
mean = plogis(beta_mu_intercept)),
aes(color = no_quota_title)) +
stat_function(fun = dprop, size = 1,
args = list(size = beta_phi,
mean = plogis(beta_mu_intercept + beta_mu_quota)),
aes(color = quota_title)) +
scale_fill_viridis_d(option = "plasma", end = 0.8, name = "Quota",
guide = guide_legend(ncol = 1, order = 1)) +
scale_color_viridis_d(option = "plasma", end = 0.8, direction = -1, name = "",
guide = guide_legend(reverse = TRUE, ncol = 1, order = 2)) +
labs(x = "Proportion of women in parliament") +
theme_clean() +
theme(legend.position = "bottom")
marginaleffects(model_beta,
newdata = typical(quota = c(FALSE, TRUE)))
augment(model_beta, newdata = tibble(quota = c(FALSE, TRUE)))
model_beta_bayes <- brm(
bf(prop_fem ~ quota,
phi ~ quota),
data = vdem_2015_fake0,
family = Beta(),
chains = 4, iter = 2000, warmup = 1000,
cores = 4, seed = 1234,
# Use the cmdstanr backend for Stan because it's faster and more modern than
# the default rstan You need to install the cmdstanr package first
# (https://mc-stan.org/cmdstanr/) and then run cmdstanr::install_cmdstan() to
# install cmdstan on your computer.
backend = "cmdstanr"
)
?set_cmdstan_path
install.packages('cmdstanr')
install.packages("cmdstanr")
model_beta_bayes <- brm(
bf(prop_fem ~ quota,
phi ~ quota),
data = vdem_2015_fake0,
family = Beta(),
chains = 4, iter = 2000, warmup = 1000,
cores = 4, seed = 1234,
# Use the cmdstanr backend for Stan because it's faster and more modern than
# the default rstan You need to install the cmdstanr package first
# (https://mc-stan.org/cmdstanr/) and then run cmdstanr::install_cmdstan() to
# install cmdstan on your computer.
backend = "cmdstanr"
)
# install.packages("remotes")
remotes::install_github("stan-dev/cmdstanr")
model_beta_bayes <- brm(
bf(prop_fem ~ quota,
phi ~ quota),
data = vdem_2015_fake0,
family = Beta(),
chains = 4, iter = 2000, warmup = 1000,
cores = 4, seed = 1234,
# Use the cmdstanr backend for Stan because it's faster and more modern than
# the default rstan You need to install the cmdstanr package first
# (https://mc-stan.org/cmdstanr/) and then run cmdstanr::install_cmdstan() to
# install cmdstan on your computer.
backend = "cmdstanr"
)
##Pacotes necessários
packages <- c("dplyr","faux")
if (any(installed_packages == FALSE)) {
install.packages(packages[!installed_packages])
}
#Carregando os pacotes
invisible(lapply(packages, library, character.only = TRUE))
install.packages("faux")
# GRAFICO 1 - Simulando a media de altura dos brasileiros
# Define os parametros usados (podem ser alterados)
renda_media <- 2700
desvio <- 2200
sims <- 100000
amostra <- 1000
renda_mediana <- 1500
# Roda a simulacao
resultado <- replicate(sims, mean(rnorm(amostra, renda_media, renda_mediana desvio)))
# Roda a simulacao
resultado <- replicate(sims, mean(rnorm(amostra, renda_media, renda_mediana, desvio)))
# GRAFICO 1 - Simulando a media de altura dos brasileiros
# Define os parametros usados (podem ser alterados)
renda_media <- 2700
desvio <- 2200
sims <- 100000
amostra <- 1000
renda_mediana <- 1500
# Roda a simulacao
resultado <- replicate(sims, mean(rnorm(amostra, renda_media, renda_mediana, desvio)))
# Roda a simulacao
resultado <- replicate(sims, mean(rnorm(amostra, renda_media, desvio)))
qplot(resultado, geom = "histogram", binwidth = 1, ylab = "Frequência", xlab = "Média de renda",
colour = I("gray"), fill = I("gray90")) +
geom_vline(xintercept = mean(resultado), size = 0.7, linetype = 2) +
scale_y_continuous(expand = c(0, 0))
# Carrega os pacotes necessarios
library(stargazer)
library(gridExtra)
library(ggplot2)
library(MASS)
qplot(resultado, geom = "histogram", binwidth = 1, ylab = "Frequência", xlab = "Média de renda",
colour = I("gray"), fill = I("gray90")) +
geom_vline(xintercept = mean(resultado), size = 0.7, linetype = 2) +
scale_y_continuous(expand = c(0, 0))
renda_media <- 2700
desvio <- 3000
sims <- 100
amostra <- 10000
# Roda a simulacao
resultado <- replicate(sims, mean(rnorm(amostra, renda_media, desvio)))
qplot(resultado, geom = "histogram", binwidth = 1, ylab = "Frequência", xlab = "Média de renda",
colour = I("gray"), fill = I("gray90")) +
geom_vline(xintercept = mean(resultado), size = 0.7, linetype = 2) +
scale_y_continuous(expand = c(0, 0))
renda_media <- 2700
desvio <- 3000
sims <- 1000
amostra <- 10000
# Roda a simulacao
resultado <- replicate(sims, mean(rnorm(amostra, renda_media, desvio)))
qplot(resultado, geom = "histogram", binwidth = 1, ylab = "Frequência", xlab = "Média de renda",
colour = I("gray"), fill = I("gray90")) +
geom_vline(xintercept = mean(resultado), size = 0.7, linetype = 2) +
scale_y_continuous(expand = c(0, 0))
renda_media <- 2700
desvio <- 40000
sims <- 1000
amostra <- 10000
# Roda a simulacao
resultado <- replicate(sims, mean(rnorm(amostra, renda_media, desvio)))
qplot(resultado, geom = "histogram", binwidth = 1, ylab = "Frequência", xlab = "Média de renda",
colour = I("gray"), fill = I("gray90")) +
geom_vline(xintercept = mean(resultado), size = 0.7, linetype = 2) +
scale_y_continuous(expand = c(0, 0))
# Roda a simulacao
resultado <- replicate(sims, median(rnorm(amostra, renda_media, desvio)))
qplot(resultado, geom = "histogram", binwidth = 1, ylab = "Frequência", xlab = "Média de renda",
colour = I("gray"), fill = I("gray90")) +
geom_vline(xintercept = mean(resultado), size = 0.7, linetype = 2) +
scale_y_continuous(expand = c(0, 0))
renda_media <- 2700
desvio <- 4000
sims <- 1000
amostra <- 10000
# Roda a simulacao
resultado <- replicate(sims, median(rnorm(amostra, renda_media, desvio)))
qplot(resultado, geom = "histogram", binwidth = 1, ylab = "Frequência", xlab = "Média de renda",
colour = I("gray"), fill = I("gray90")) +
geom_vline(xintercept = mean(resultado), size = 0.7, linetype = 2) +
scale_y_continuous(expand = c(0, 0))
summary(resultado)
renda_media <- 2800
desvio <- 4000
sims <- 1000
amostra <- 10000
# Roda a simulacao
resultado <- replicate(sims, median(rnorm(amostra, renda_media, desvio)))
summary(resultado)
# Roda a simulacao
resultado <- replicate(sims, median(rbinom(amostra, renda_media, desvio)))
summary(resultado)
# Roda a simulacao
resultado <- replicate(sims, median(rnorm(amostra, renda_media, desvio)))
summary(resultado)
qplot(resultado, geom = "histogram", binwidth = 1, ylab = "Frequência", xlab = "Média de renda",
colour = I("gray"), fill = I("gray90")) +
geom_vline(xintercept = mean(resultado), size = 0.7, linetype = 2) +
scale_y_continuous(expand = c(0, 0))
summary(resultado)
rand_unifs_10 <- runif(n = 10, min = 0, max = 1)
hist(rand_unifs_10000, xlab = "Random value (X)", col = "grey",
main = "", cex.lab = 1.5, cex.axis = 1.5)
rand_unifs_10000 <- runif(n = 10000, min = 0, max = 1);
rand_unifs_10000 <- runif(n = 10000, min = 0, max = 1)
hist(rand_unifs_10000, xlab = "Random value (X)", col = "grey",
main = "", cex.lab = 1.5, cex.axis = 1.5)
rand_poissons <- rpois(n = 10, lambda = 1.5)
print(rand_poissons)
rand_poissons_10000 <- rpois(n = 10000, lambda = 4.5)
hist(rand_poissons_10000, xlab = "Random value (X)", col = "grey",
main = "", cex.lab = 1.5, cex.axis = 1.5)
mean(rand_poissons_10000)
median(rand_poissons_10000)
rand_poissons_10000_2 <- rpois(n = 100000, lambda = 7.5)
hist(rand_poissons_10000_2, xlab = "Random value (X)", col = "grey",
main = "", cex.lab = 1.5, cex.axis = 1.5)
mean(rand_poissons_10000_2)
median(rand_poissons_10000_2)
vecpoisson=rpois(100,5)
mean(vecpoisson)
set.seed(198911)
vecpoisson=rpois(100,5)
mean(vecpoisson)
reps <- 50000
nexps <- 5
rate <- 0.1
set.seed(0)
system.time(
x1 <- replicate(reps, sum(rexp(n=nexps, rate=rate)))
) # replicate
summary(x1)
require(ggplot2)
ggplot(data.frame(x1), aes(x1)) +
geom_histogram(aes(y=..density..)) +
stat_function(fun=function(x)dgamma(x, shape=nexps, scale=1/rate),
color="red", size=2)
qplot(resultado, geom = "histogram", binwidth = 1, ylab = "Frequência", xlab = "Média de renda",
colour = I("gray"), fill = I("gray90")) +
geom_vline(xintercept = mean(resultado), size = 0.7, linetype = 2) +
scale_y_continuous(expand = c(0, 0))
setwd("C:/Users/Bruno Schaefer/Documents/IESP/Editoria Dados/Papers/João What is liberalism/data/2 Quant data analysis")
# Carregar o pacote "foreign"
library(foreign)
# Ler o arquivo Stata
dados <- read.dta("PELA_LAPOP.dta")
# Ler o arquivo Stata
dados <- read.dta("PELA_LAPOP.dta")
# Carregar o pacote "foreign"
library(haven)
# Ler o arquivo Stata
dados <- read_dta("PELA_LAPOP.dta")
View(dados)
View(dados)
summary(dados)
# Criar a tabela de contingência
tabela5 <- table(dados$parties, dados$country)
# Calcular os percentuais nas células da tabela
tabela_percentual5 <- prop.table(tabela5) * 100
tabela_percentual5
# Criar um data frame com as variáveis
tabela5 <- data.frame(dados$parties, dados$country, dados$publicfirms)
# Calcular a média da variável var3 por combinação das variáveis var1 e var2
media5 <- aggregate(publicfirms ~ parties + country, data = tabela5, FUN = mean)
summary(dados$publicfirms)
View(tabela5)
# Calcular a média da variável var3 por combinação das variáveis var1 e var2
media5 <- aggregate(publicfirms ~ parties + country, data = dados, FUN = mean)
# Exibir a tabela com a média da variável var3
print(media5)
# Instalar e carregar os pacotes necessários
install.packages(c("knitr", "kableExtra"))
library(knitr)
library(kableExtra)
install.packages(c("knitr", "kableExtra"))
# Exibir a tabela com a média da variável var3
media5 <- media5 %>% kable()
kableExtra::kable_styling(media5)
table(dados$database)
dados_eleitores <- subset(dados, c(dados$database == 2))
dados_manifestos <- subset(dados, c(dados$database == 1))
# Criar a tabela de contingência
tabela_n <- table(dados_manifestos$parties, dados_manifestos$country)
# Calcular os percentuais nas células da tabela
tabela_percentualn <- prop.table(tabela_n) * 100
tabela_percentualn
summary(dados_manifestos$publicfirms)
# Calcular a média da variável var3 por combinação das variáveis var1 e var2
media_n <- aggregate(publicfirms ~ parties + country, data = dados_manifestos, FUN = mean)
# Exibir a tabela com a média da variável var3
print(media_n)
# Exibir a tabela com a média da variável var3
media_n <- media_n %>% kable()
kableExtra::kable_styling(media_n)
# Criar a tabela de contingência
tabela_s <- table(dados_eleitores$parties, dados_eleitores$country)
# Calcular os percentuais nas células da tabela
tabela_percentuals <- prop.table(tabela_s) * 100
tabela_percentuals
summary(dados_eleitores$publicfirms)
# Calcular a média da variável var3 por combinação das variáveis var1 e var2
media_s <- aggregate(publicfirms ~ parties + country, data = dados_eleitores, FUN = mean)
# Exibir a tabela com a média da variável var3
print(media_s)
# Exibir a tabela com a média da variável var3
media_s <- media_s %>% kable()
kableExtra::kable_styling(media_s)
summary(dados_manifestos$democracy)
# Calcular a média da variável var3 por combinação das variáveis var1 e var2
media_x <- aggregate(democracy ~ parties + country, data = dados_manifestos, FUN = mean)
# Exibir a tabela com a média da variável var3
print(media_x)
# Exibir a tabela com a média da variável var3
media_x <- media_x %>% kable()
kableExtra::kable_styling(media_x)
summary(dados_eleitores$democracy)
summary(dados_eleitores$democracy)
summary(dados_eleitores$democracy_cont)
# Calcular a média da variável var3 por combinação das variáveis var1 e var2
media_z <- aggregate(democracy_cont ~ parties + country, data = dados_eleitores, FUN = mean)
# Exibir a tabela com a média da variável var3
print(media_z)
# Exibir a tabela com a média da variável var3
media_z <- media_z %>% kable()
kableExtra::kable_styling(media_z)
View(dados_eleitores)
# Exibir a tabela com a média da variável var3
print(media_z)
# Exibir a tabela com a média da variável var3
print(media_z)
# Exibir a tabela com a média da variável var3
media_z <- media_z %>% kable()
kableExtra::kable_styling(media_z)
summary(dados_eleitores$democracy)
summary(dados_eleitores$democracy_cont)
View(dados_eleitores)
summary(dados_manifestos$samesex)
# Calcular a média da variável var3 por combinação das variáveis var1 e var2
media_r <- aggregate(samesex ~ parties + country, data = dados_manifestos, FUN = mean)
# Exibir a tabela com a média da variável var3
print(media_r)
# Exibir a tabela com a média da variável var3
media_r <- media_r %>% kable()
kableExtra::kable_styling(media_r)
summary(dados_eleitores$samesex)
# Calcular a média da variável var3 por combinação das variáveis var1 e var2
media_ss <- aggregate(samesex ~ parties + country, data = dados_eleitores, FUN = mean)
# Exibir a tabela com a média da variável var3
print(media_ss)
# Exibir a tabela com a média da variável var3
media_ss <- media_ss %>% kable()
kableExtra::kable_styling(media_ss)
summary(dados_manifestos$reduceinequality)
# Calcular a média da variável var3 por combinação das variáveis var1 e var2
media_ri <- aggregate(reduceinequality ~ parties + country, data = dados_manifestos, FUN = mean)
# Exibir a tabela com a média da variável var3
print(media_ri)
# Exibir a tabela com a média da variável var3
media_ri <- media_ri %>% kable()
kableExtra::kable_styling(media_ri)
summary(dados_eleitores$reduceinequality)
# Calcular a média da variável var3 por combinação das variáveis var1 e var2
media_rin <- aggregate(reduceinequality ~ parties + country, data = dados_eleitores, FUN = mean)
# Exibir a tabela com a média da variável var3
print(media_rin)
# Exibir a tabela com a média da variável var3
media_rin <- media_rin %>% kable()
kableExtra::kable_styling(media_rin)
summary(dados_eleitores$democracy)
summary(dados_eleitores$democracy_cont)
# Calcular a média da variável var3 por combinação das variáveis var1 e var2
media_z <- aggregate(democracy_cont ~ parties + country, data = dados_eleitores, FUN = mean)
# Exibir a tabela com a média da variável var3
print(media_z)
# Exibir a tabela com a média da variável var3
media_z <- media_z %>% kable()
kableExtra::kable_styling(media_z)
table(dados_manifestos$parties, dados_manifestos$country)
table(media_n)
# Exibir a tabela com a média da variável var3
print(media_n)
# Calcular a média da variável var3 por combinação das variáveis var1 e var2
media_n <- aggregate(publicfirms ~ parties + country, data = dados_manifestos, FUN = mean)
table(media_n)
# Exibir a tabela com a média da variável var3
print(media_n)
